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SUMMARY 


A  new  and  simple  method  of  finite-element  grid  improvement  is  presented. 
The  objective  is  to  improve  the  accuracy  of  the  analysis.  The  procedure  is 
based  on  a  minimization  of  the  trace  of  the  stiffness  matrix.  For  a  broad 
class  of  problems  this  minimization  is  seen  to  be  equivalent  to  minimizing  the 
potential  energy.  The  method  is  illustrated  with  the  classical  tapered  bar 
problem  examined  earlier  by  Prager  and  by  Masur.  Identical  results  are 
obtained . 


INTRODUCTION 

In  a  general  context,  the  finite-element  method  is  an  approximate  proce¬ 
dure  for  solving  differential  equations.  The  accuracy  of  the  method  depends 
on  (1)  the  number  of  elements,  (2)  the  choice  of  interpolation  functions,  and 
(3)  the  location  of  the  grid  points.  The  number  of  elements,  and  hence  the 
number  of  grid  points,  is  usually  restricted  simply  by  computer  capability  and 
by  processing  costs.  Also,  there  have  been  many  recent  advances  in  improving 
the  accuracy  of  the  finite-element  method  by  using  higher  order  interpolation 
polynomials  and  shape  functions  (p-method)  and  by  exhaustive  analysis  with 
large  numbers  of  elements  (h-method).  Some  of  these  advances  are  described  in 
the  references  cited  herein.  However,  optimal  grid  point  location  (r- method) 
is  far  less  advanced.  Practical  procedures  for  the  analyst  still  need  to  be 
developed  and  refined. 

One  of  the  earliest  attempts  to  develop  a  grid  optimization  procedure  was 
that  of  Prager  in  1975  (ref.  1).  Prager's  work  provided  a  stimulus  and  a  basis 
for  later  grid  optimization  research  as  recorded  in  references  2  to  19.  Note¬ 
worthy  among  these  efforts  are  the  works  of  Shepard  (refs.  11,  13,  and  15), 
Masur  (ref.  2),  Turcke  (refs.  8  and  9),  Carroll  (ref.  7),  NcNiece  (ref.  4), 
Carey  (ref.  16),  Diaz  (ref.  12),  Melosh  (ref.  10),  Durocher  (ref.  17),  and 
their  col  leagues . 


Prager  examined  a  bar  with  a  linearly  varying  cross  section  under  tension. 
He  showed  that  the  grid  producing  the  desired  least  potential  energy  is  the  one 
where  the  cross-section  areas  at  the  nodes  form  a  geometric  series.  In  this 
configuration,  the  strain  energy  is  divided  equally  among  the  elements. 
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Masur  (ref.  2)  observes  that  this  latter  result  of  equal  element  strain 
energies  is  not  a  general  characteristic  of  optimal  meshes  but  instead  is  a 
result  of  the  simple  geometry  of  Prager's  problem. 

In  this  paper  we  present  a  finite-element  grid  improvement  technique  which 
is  based  on  the  minimization  of  the  trace  of  the  global  stiffness  matrix.  We 
show  that  this  method  leads  to  identical  results  to  those  of  Prager.  It  has 
the  advantage  of  being  simpler  than  traditional  optimization  procedures. 

The  method  presented  herein  provides  a  mesh  improvement  which  is  based  on 
the  geometry  of  the  body.  As  such,  it  provides  a  significant  improvement  over 
uniform  meshes,  and  it  produces  a  good  first  iteration  for  accommodating  spe¬ 
cial  loading  configurations. 

In  the  usual  finite-element  procedure,  the  governing  equations  are 
obtained  by  minimizing  a  functional  tt  by  varying  the  dependent  variables  of 
the  physical  problem  (ref.  20).  For  el astostati cs  this  is  equivalent  to  the 
principle  of  minimum  potential  energy  (ref.  21).  This  leads  to  the  familiar 
system  of  linear  algebraic  equations.  Attempts  to  minimize  ir  with  respect 
to  the  nodal  coordinates,  however,  leads  to  a  system  of  nonlinear  equations. 
These  equations  are  generally  extremely  difficult  to  solve  even  for  the  sim¬ 
plest  cases.  To  avoid  this  difficulty  we  are  proposing  instead  to  examine  the 
stiffness  matrix  to  obtain  information  about  an  improvement  in  grid  point  loca¬ 
tion.  Our  motivation  is  the  observation  of  the  major  role  of  the  stiffness 
matrix  in  the  value  of  the  potential  ir.  Also,  the  entries  of  the  stiffness 
matrix  are  dependent  on  the  grid  point  coordinates. 


NOMENCLATURE 


A^  nodal  area 

A g_  end  area 

A0  base  area 

A|<  area  of  element 

c  area  ratio  (see  eq.  (7)) 

E  elastic  modul us 

{f}  global  force  array 

{f }  transformed  global  force  array 

fi  entries  of  {f } 

CK]  global  stiffness  matrix 
C K ]  di agonal  form  of  [K] 


nCvCv-:-/vCv->,cvc 
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k.  entries  of  [K] 

L  bar  length 

Q_k  length  of  k^  element 

n  number  of  elements 

P  axial  load 

r i  j  length  ratio, 

n 

S  sum  of  length  ratios,  £  r.. 

"  3-1  Jl 

T  trace  of  C K ] 

[T]  transformation  matrix 

u  axial  displacement 

{u}  global  displacement  array 
{u}  transformed  global  displacement  array 

Uj  entries  of  {u } 

x  axial  coordinate 

x^  nodal  coordinate 

k  element  stiffness  matrix 

it  potential  energy 

5  dimensionless  axial  coordinate 

5k  dimensionless  nodal  coordinate 

ANALYSIS 

Our  objective  is  to  develop  a  practical  and  efficient  procedure  of  grid 
enhancement  tending  towards  optimization.  Our  thesis  is  that  for  many  problems 
the  minimization  of  the  trace  of  the  stiffness  matrix  leads  to  a  minimization 
of  the  potential  energy  and,  as  a  consequence,  provides  the  optimal  grid  con- 
f i gurat i on . 

To  see  this,  consider  the  governing  matrix  equation  of  finite-element 
analysi s : 

[ K ]  { u >  =  {f }  (1) 

where  [K]  is  the  stiffness  matrix,  {u}  the  array  of  dependent  variables,  and 
{f}  the  force  array.  We  can  view  [K]  as  an  operator  which  maps  { u }  into 


(f}.  In  this  context,  since  [K]  is  symmetric1  we  can  find  an  orthogonal 
transformation  [T]  which  diagonalizes  [ K ] ;  that  is. 


CK]  =  [T]T[K][T] 


(2) 


where  [K]  is  a  diagonal  matrix.  Let  C T 3 { u }  and  [ T ] { f >  be  {u}  and  { f } . 
Then  the  potential  energy  ir  may  be  expressed  as 


Tt  =  \  { u } T  [  K  ]  {u}  -  {f  }T{U}  =  |  { u } T  C  K  ]  {u}  -  {f}T{u}  (3) 


In  terms  of  the  array  components,  ir  becomes 


IT  = 


n 


i  =  1 


(4) 


where  the  kj  (i  =  l,...,n)  are  the  diagonal  entries  of  CK] . 

n  „  ^ 

Observe  in  equation  (4)  that  the  last  term  -  E  f.u.  does  not  explicitly 

U1  1  1 
n  ^  „ 

involve  the  nodal  coordinates.  Therefore,  -  E  f.u.  does  not  affect  the 

i  =  1 

minimization  of  tt  with  respect  to  the  nodal  coordinates.  Also,  since  the  u] 
are  positive  and  are  independent  variables  in  the  minimization  of  ir,  the  mini¬ 
mization  of  ir  with  respect  to  the  nodal  coordinates  occurs  when  the  sum  of 

the  kj  (the  trace  of  [K])  is  a  minimum.  Since  the  trace  of  a  matrix  is 

invariant  under  an  orthogonal  transformation,  minimizing  the  trace  of  C K ]  is 
equivalent  to  minimizing  the  trace  of  [XL 


In  minimizing  the  trace,  we  will  not  adversely  affect  the  diagonal  domi¬ 
nance  of  [K]  required  to  avoid  i 1 1 -condi t ion i ng .  The  improved  stiffness 
matrix  we  seek  is  the  result  of  redistribution  of  the  nodes  and  not  of  an  arbi¬ 
trary  mathematical  operation. 

To  illustrate  the  application  of  these  concepts,  consider  the  axially 
loaded  tapered  bar  shown  in  figure  1.  (This  is  the  same  problem  examined  by 
Prager  (ref.  1)  and  Masur  (ref.  2).)  The  objective  is  to  determine  a  finite- 
element  mesh  which  best  predicts  the  axial  displacement.  Let  the  bar  have 
length  L  and  let  it  be  divided  into  n  elements  with  n  +  1  nodes  (numbered 
0  to  n)  as  shown.  Let  the  areas  at  the  ends  of  the  bar  be  A0  and  A^.  Let 
f;  be  the  nond i mens i ona  1  length  parameter  defined  as 


(5) 


’The  analysis  which  follows  is  based  on  the  symmetry  of  [Kj.  If  [K] 
is  not  symmetric,  a  similar  analysis  could  be  developed  using  nonorthogonal 
transformat  ions. 


Then  the  area  at  any  particular  E,  along  the  bar  is 


A  =  A0( 1  -  c£) 


where  c  is 


A  -  An 
o  1 


o  <  c  <  1 


Hence,  the  area  at  the  kth  node  is 


A|<  =  A0(  1  -  cEk ^ 


where  Ek  is  £(xk>. 


Let  the  individual  elements  have  a  uni form_cross  section.  For  example, 
let  the  k^h  element  have  cross-section  area  Ak  and  length  &k  as  in  fig¬ 
ure  2.  (Note  that  the  elements  do  not  necessarily  have  the  same  length.)  Then 
Ak  and  2k  are 


A  +  A 
k-1  k 


2k  =  *k  "  xk-l  =  L(^k  "  ^k-l} 


The  element  stiffness  matrix  for  the  k^h  element  is  (ref.  20) 


A.E  1  -1 

[,1  =  ±- 

k  -1  1 


where  E  is  the  elastic  modulus.  Then  the  trace  T  of  the  global  stiffness 
matrix  is 


T  =  2E 


_1  E 

2  i  L 


(k.  ,  +  A. 

l«t  -  *1-1 


The  improved  node  location  occurs  when  the  trace  is  minimized  with  respect 
to  the  nodal  coordinates  Ek(k=l , .  . .  ,n-l )  .2  Hence,  by  setting  the  derivative 
of  T  with  respect  to  Ek  equal  to  zero,  we  obtain 

^Ncdes  numbers  0  and  n  are  fixed  at  the  ends  of  the  bar  and  thus  not 
considered  for  optimization. 


(A,  +  A 


A.  +  A 


Using  equation  (9)  and  simplifying  we  obtain 


A  i  ,  =  A, 
k+ 1  k 


k-1  a 


To  simplify  the  analysis  it  is  convenient  to  introduce  the  length  ratio 
parameter  r^j  defined  as  2j/ftj.  Then  the  ratio  ftk+l/Q-k  may  be  written  as 


Then  L/ft.  is 


k+1  ,1 


rkl  rkl 


where  S  is  defined  as  E  r...  Using  this  notation,  we  can  rewrite 
n  ,  i  1 

equation  (13)  as 


A,  .  =  A,  P±Li 

k+1  k  l  rkl 


k-1  r 


rk+l , A  cAork+l , 1  (r k+1,1 


Also,  L  may  be  written  as 


2.  +  ft,  + 


1  +  r  „  . 


r 


i  ki'L*  h*  V  k*  It*1  A#  VUi1  .V  .M.fc*  tft'.V  Jt*  J  .ViV  AfAVt* 


v.v^.v.v; 


Then,  from  equation  (8),  A,  may  be  written  as 


A.  =  A 
k  o 


Specifically,  A^  and  A^  are 


A.  =  A  1  - 
1  o 


A2  -  Ao  1  -  c 


(I  .  r 2 1 ) 


To  obtain  the  element  area  ratios  let  k  =  1  in  equation  (16).  A^  is  then 

f  2  \  2  r  2 1  +  ^ 

A2  =  Al^r21"  ]J  +  Aor21  +  cAor21  ~T^  (' 

Then,  by  using  equations  (19)  to  eliminate  A]  and  we  obtain 

Sn( 1  "  r21 }  =  c 

From  the  first  of  equations  (19)  we  have 

A,  =  A  r_ . 

1  o  21 


A  ”  r21 
o 


Similarly,  the  second  of  equations  (19)  leads  to 


A-  =  A  rt. , 
2  o  21 


A,  =  '21 


V  N  N  V  S  ^.  V*  ^  ^  ^ 


Next,  let  k  =  2  in  equation  (16).  By  using  the  same  procedure  we  obtain 


r32  "  S  (1  r32  +  r21} 

n 


But,  in  view  of  equation  (21)  this  becomes 


or 


Thus , 


1  r32  -  (1  r2l )(1  r 32  +  r21 ) 


r21(r3 2  r21)  =  0 


r32  "  r2 1 


From  equation  (18)  we  have 


J 


A,.  =  A 
3  o 


fh' 


=  A  r 


o  21 


Therefore,  we  obtain 


A3  A2  A^ 

AT  =  A?  =  A-  =  r21  =  r32 

2  1  o 


Proceeding  similarly  for  k  =  3,4,..,  we  obtain 

Vi  ^ 

Vi =  v2  =  • =  A0  =  "21  =  r32  = ' 


=  r 


n ,  n-1 


Thus,  we  have  the  relations 


r31  =  r  2 1 r  32  =  r21 


41  "  43  3221  "  21 


k-1 
kl  ‘  '21 


(24) 


(25) 


(26) 


(27) 


(28) 


(29) 


(30) 


Hence,  S,  is  the  geometric  series 


Then,  from  equations  (18)  and  (21),  we  have 


A.  -  A 
k  o 


1  -  (1  -  r,,>  =  A  r 


21 


)  = 


O  21 


and  then 


A  =  A  r».  =  An 
n  o  21  fi. 


(32) 


Then,  from  equation  (7)  we  see  that  r^  is 


r21  =  ^ "  c) 


1/n 


(33) 


Finally,  by  substituting  into  equation  (17)  we  have 


k-1 


Ek  “  L  (l  +  r2 1  +  r31  + 


+  rkl}  =  (1  -  r21} 


1  +  r21  +  r21  +  '  '  '  +  r21 


1  "  r21  1  -  (1  -  c)k/n 


(34) 


This  is  the  result  obtained  by  Prager  (ref.  1)  in  his  analysis  of  the  same 
problem. 


DISCUSSION 

First,  observe  that  in  equation  (34)  for  a  uniform  thickness  beam  c  =  0 
and  thus  Fk  's  undetermined.  This  means  that  for  a  uniform  thickness  beam 
the  nodal  positions  are  arbitrary;  that  is,  all  meshes  are  equally  optimal  for 
a  uniform  thickness  beam. 

Next,  consider  again  the  element  stiffness  matrix  of  equation  (10).  From 
equations  (8)  and  (34)  the  scalar  multiplier  is 


EAk  E(Ak-i  +  V  Aoc  +  r2A 

ftk  =  2L(Ek  -  Ek-1>  =  2L  l1  -  r2 1  / 


(35) 


Since  this  is  a  constant  (independent  of  k)  the  element  stiffness  matrix  is 
the  same  *cr  each  element.  This  means  that  each  element  has  the  same  strain 
energy.  Masur  (ref.  2)  has  suggested  that  this  result  is  due  to  the  simple 
■geometry  :f  the  problem. 

Even  with  this  simple  geometry,  however,  the  analysis  needed  to  determine 
the  optimal  nodal  positions  has  been  extremely  detailed.  With  more  complex 
geometries  the  analysis  will  become  intractable.  Alternatively,  a  more  conven¬ 
ient  method  of  improving  the  nodal  positions  is  to  examine  the  trace  of  the 
stiffness  matrix  and  to  adjust  the  nodal  positions  to  minimize  the  trace.  The 
criteria  for  minimizing  the  trace  of  the  stiffness  matrix  is  a  comparatively 
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simple  procedure  -  readily  amenable  to  the  development  of  computer  algorithms 
for  optimal  nodal  locations. 


In  the  preceeding  discussion,  we  have  ignored  any  cons ideration  of  the 
effects  of  concentrated  loads  or  boundary/edge  effects.  It  is  common  engineer¬ 
ing  practice  to  use  a  fine  (dense)  grid  in  highly  loaded  regions  and  to  place 
a  grid  point  specifically  at  the  point  of  application  of  a  concentrated  load. 
The  trace  minimization  suggested  here  is  intended  to  aid  in  the  discretization 
process  at  locations  removed  from  concentrated  loads.  Our  justification  is 
based  on  Saint-Venant ' s  principle  that  localized  effects  disappear  at  short 
di stances . 

A  principal  question  with  the  minimum  trace  method  is:  What  is  the  range 
of  applicability?  The  authors  have  found  the  method  to  lead  to  significant 
decreases  in  potential  energy  from  that  of  a  uniform  mesh  for  structural  and 
heat-transfer  problems.  The  range  of  applicability  is  currently  being 
explored.  Numerical  algorithms  using  this  procedure  are  being  developed. 
Finally,  the  influence  of  values  of  second  and  higher  invariants  of  the  stiff¬ 
ness  matrix  needs  to  be  explored. 


NUMERICAL  EXAMPLE 


To  illustrate  the  value  of  optimizing  the  mesh,  consider  an  axially  loaded 
bar  which  tapers  to  1/3  the  base  area  as  in  figure  3.  Specifically,  let  P, 

A0,  c,  E,  and  L  have  the  following  values: 

P  =  20  N  'l 

Aq  =  0.0015  m2 


E  =  2 . 07x 1 01 1  N/m2 
L  =  4  m 

The  objective  is  to  find  the  axial  displacement. 


(35) 


J 


From  elementary  mechanics  the  axial  displacement  u  at  any  location  < 
i  s 


u 


PL 

A  Ec 
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r 


<  37) 


To  compa  e  *:he  displacement  results  of  finite-element  mode’s  with  equa¬ 
tion  (37),  four  models  of  the  bar,  each  having  four  elements,  we,-e  examined. 

One  :f  the  models  had  a  uniform  nodal  distribution.  Another  had  the  'optimal" 
me sn  as  developed  in  equation  (34).  The  remaining  two  models  had  arbitarily 
se’ectea  nodal  distributions.  The  nodal  displacements  were  evaluated  using  tne 
mode’s  and  compa;  ed  with  the  displacement  calculated  by  equation  (3?’1. 
Tab’e  I  stows  the  results.  T a b l e  II  presents  an  error  analysis  and  also  an 
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fulness  of  the  trace  minimization  mesh  improvement  method. 

2.  Minimization  of  the  trace  of  the  stiffness  matrix  is  a  relative';/ 
simple  mesh  optimization  procedure.  It  is  readily  adaptable  to  algorithm 
development. 

3.  Trace  minimization  can  be  used  in  combination  with  other  grid  optimiza¬ 
tion  tecnniques.  Indeed,  it  can  provide  a  starting  mesh  for  iteration  tech¬ 
niques,  useful  for  specialized  loading. 

4.  The  principal  benefit  of  the  trace  minimization  method  is  accuracy  of 
analysis  as  opposed  to  efficiency  of  analysis. 
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TABLE  I.  -  COMPARISON  OF  AXIAL  DISPLACEMENTS  FOR  THE  BAR  OF  FIGURE  3  CALCULATED 

USING  VARIOUS  MODELS 


Axial 
location, 
x » 

m 

Exact  displacement 
(eq.  (37)). 

10"9  m 

Displacements  computed^usi  ng  various  models, 
10-9  m 

Uni  form 
mesh 

"Optimum" 

mesh 

(Prager) 

Mesh  3 

Mesh  4 

0.0 

5 

0.0 

33.6276 

70.46244 

106.1461 

156.702 

208.3078 

212.2923 

267.883 

318.4384 

424.5845 

0.0 

0.0 

0.0 

33.60639 

70.41338 

0.0 

33.60639 

ilo 

1.4409860 

2.0 

2.5 

2.5358986 

3.0 

3.3678522 

4.0 

70.2679 

105.482 

156.1509 

15.56506 

207.1804 

206.8158 

210.9648 

266.5719 

316.4529 

421.940 

421 . 1612 

417.6195 

417.9814 

TABLE  II.  -  ERROR  ANALYSIS 


Axial 
location, 
x « 

m 

Error  for  various  meshes, 

10“9  m 

Uni  form 
mesh 

"Optimum" 

mesh 

(Prager) 

Mesh  3 

Mesh  4 

0.0 

.5 

1.0 

1 .441 

2.0 

2.5 

2.536 

3.0 

3.368 

4.0 

L2~Norm 

0.0 

0.0 

.02121 

.04906 

0.0 

.02121 

.1945411 

.664118 

.5506105 

1.0514 

1.1274 

1.492 

1.32769 

1.311108 

1 .985439 

2.64433 

3.6246009 

3.423228 

3.7119418 

6.965 

7.1232118 

6.6031 

6.780697 

FIGURE  1.  -  LINEAR  TAPERED  BAR. 
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